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ABSTRACT 

We use the NIHAO (Numerical Investigation of Hundred Astrophysical Objects) cos¬ 
mological simulations to investigate the effects of baryonic physics on the time evolution 
of Dark Matter central density profiles. The sample is made of Ri 70 independent high 
resolution hydrodynamical simulations of galaxy formation and covers a wide mass range: 
10^° < Mhaio/M 0 ;< 10^^, i.e., from dwarfs to L*. We confirm previous results on the 
dependence of the inner dark matter density slope, a, on the ratio between stellar-to-halo 
mass, Mstar/Afhaio- We show that this relation holds approximately at all redshifts (with an 
intrinsic scatter of ^ 0.18 in a measured between 1 — 2% of the virial radius). This implies 
that in practically all haloes the shape of their inner density prohle changes quite substan¬ 
tially over cosmic time, as they grow in stellar and total mass. Thus, depending on their final 
Afstar/Afhaio ratio, haloes can either form and keep a substantial density core (i?core 1 
kpc), or form and then destroy the core and re-contract the halo, going back to a cuspy profile, 
which is even steeper than CDM predictions for massive galaxies ( 1 O^^M 0 ). We show that 
results from the NIHAO suite are in good agreement with recent observational measurements 
of a in dwarf galaxies. Overall our results suggest that the notion of a universal density profile 
for dark matter haloes is no longer valid in the presence of galaxy formation. 
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1 INTRODUCTION 

The Cold Dark Matter (CDM) theory is very successful at de¬ 
scribing the Universe’s topology and evolution on large scales (e.g. 
Springel et al. 2005). This theory makes clear predictions of the dis¬ 
tribution of Dark Matter on small scales, where all collapsed struc¬ 
tures (haloes) are supposed to share an approximately self-similar 
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dark matter density profile, as first pointed out by Navarro, Frenk 
and White (1997, NFW) and then confirmed via higher resolution 
numerical simulations by a large variety of works (e.g., Klypin 
et al. 2001, Diemand et al. 2004; Power et al. 2003, Maccio et al 
2007, 2008, Neto et al. 2007, Navarro et al. 2010, Prada et al. 2012, 
Dutton & Maccio 2014). 

This prediction of a universal, cuspy profile for the dark mat¬ 
ter seems to be in tension with observations of rotation curves of 
low mass galaxies (e.g Moore 1994, Salucci & Burkert 2000; Si¬ 
mon et al. 2005; de Blok et al. 2008; Kuzio de Naray, McGaugh 
& de Blok 2008; Kuzio de Naray, McGaugh & Mihos 2009; Oh 
et al. 2011a), which seem to prefer density profiles with a constant 
density core in the center. 
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Table 1. Properties of our three selected galaxies at z = 0 


simulation ID 

# particles 

7^200 

# DM particles 

# Star pailicles 

^star 

DM pailicle mass 

moM 

virial mass 
(-A^2OO/M0) 

Stellar mass 

(A^star/M©) 

2 = 0 SFR 
(M© yr-l) 

Conv. rad.^ 
(kpc) 

g2.63el0 

458,723 

414,291 

18,388 

1.173 X lO'^ 

2.70 X IQi® 

4.28 X 10’^ 

0.00 

0.41 

g2.19ell 

920,447 

557,247 

113,958 

2.169 X 10® 

1.31 X IQii 

9.27 X 10* 

8.2 X 10-2 

0.65 

g8.06ell 

1,366,038 

481,349 

665,796 

1.735 X 10® 

9.43 X IQii 

4.48 X lO^® 

11.31 

1.23 


“ Following Power et al. 2003, see text for details. 


This cusp-core controversy suggests that either a modification 
of the whole CDM paradigm is required (e.g., self interacting dark 
matter, SIDM, Vogelsberger et al. 2014), or the inadequacy of pure 
N-body simulations to capture the dark matter dynamics on small 
scales, due to the absence of dissipative phenomena connected to 
baryonic physics. 

Recent and more accurate cosmological hydrodynamical sim¬ 
ulations have indeed shown that baryons are able to alter the dark 
matter profile and to create substantial cores (up to ~ 1 kpc) in 
the dark matter distribution (e.g. Governato et al. 2010, Maccio 
et al. 2012, Martizzi et al. 2013, Munshi et al. 2013, Di Cintio 
2014a,b, Trujillo-Gomez et al. 2015, Onorbe et al. 2015). 

As nicely described in Pontzen & Governato (2012) the for¬ 
mation of dark matter cores is linked to the rapid change of the 
total gravitational potential. In the central kiloparsecs the potential 
changes on sub-dynamical time-scales, as repeated, energetic feed¬ 
back generates large underdense bubbles of expanding gas. The 
fluctuations in the central potential irreversibly (e.g. non adiabat- 
ically) transfer energy into collisionless particles, thus generating a 
dark matter core (see also Ogiya & Mori 2014 for a detailed analy¬ 
sis of the effect of resonances between DM particles and the density 
wave excited by the oscillating potential). These strong gas out¬ 
flows are triggered by stellar (and AGN) feedback (Navarro, Eke 
& Frenk 1996; Read & Gilmore 2005; Mashehenko, Couchman & 
Wadsley 2006; Pontzen & Governato 2014; Martizzi et al. 2013). 

Recently Di Cintio et al. (2014a, DC14 hereafter) has shown 
that the modification of the initial dark matter profile (either lead¬ 
ing to an expansion or a contraction) is clearly linked to the in¬ 
tegrated star formation efficiency of the galaxy, which can be pa¬ 
rameterized through the redshift zero stellar mass - halo mass ra¬ 
tio (Mstar/Mhaio), and that this holds for a large variety of stellar 
feedback implementations. 

They clearly showed that the halo response to star formation 
is non monotonic and that haloes with very low and very high 
star formation efficiency (or with very low and high Mstar/Afhaio) 
tend to preserve the initial cuspy profile, while haloes with 
Afstar/Mhaio ~ 0.3% have a flat central density profile. 

In this paper we want to expand the original work of DC 14, 
by using a newer (and larger) sample of simulated galaxies. This 
new set of simulations has been performed using the new cosmo¬ 
logical parameters from the Planck satellite (the Planck Collabora¬ 
tion 2014), and an updated version of the GASOLINE code (Keller 
et al. 2014) which fixes some of the known problems of the SPH 
technique (Agertz et al. 2007). This newer and larger sample allows 
us to look at the evolution of the density profiles through cosmic 
time and to witness the creation (and destruction) of density cores 
as a function of stellar mass and halo mass. 

This paper is organized as follows: In §2 we give an overview 
of the cosmological hydrodynamical simulations used in this work. 


Results are presented in §3, including core creation in §3.2 and core 
destruction in §3.3. Our conclusions are given in §4. 


2 SIMULATIONS 

In this study we use simulations from the NIHAO suite (Wang 
et al. 2015), based on an updated version of the MaGICC project 
(Stinson et al. 2012). All the simulations have adopted the latest 
compilation of cosmological parameters from the Planck satellite 
(the Planck Collaboration 2014); namely Dm = 0.3175, Hq = 
67.1, (78 = 0.8344, n = 0.9624 and Db = 0.0490. The haloes 
to be re-simulated with baryons and higher resolution have been 
extracted from three different pure N-body simulations with a box 
size of 60, 20 and 15 /i“^Mpc respectively, more information on 
the collisionless simulations can be found in Dutton & Maccio 
(2014). 

The aim of the NIHAO project is to study galaxy formation 
over a large mass range, from dwarf galaxies to massive spirals 
such as the Milky Way. We have decided to keep the same rela¬ 
tive mass resolution accross the whole mass range, meaning that at 
high resolution we have (roughly) the same number of dark matter 
particles within the virial radius of our galaxies. This requirement 
sets the mass of the dark matter particle and the initial mass of the 
gas particles, since the latter is simply obtained by the dark matter 
mass multiplied by (Dm — Db)/Db- 

The zoomed initial conditions were created using a modified 
version of GRAFIC2 (Bertschinger 2001, Penzo et al. 2014). The 
starting redshift is Zstart = 99, and each halo is initially sim¬ 
ulated at high resolution with DM-only using PKDGRAV (Stadel 
2001). More details on the sample selection can be found in Wang 
et al. (2015). We refer to simulations with baryons as the hydro 
simulation, while we will use the term N-body for the DM-only 
simulation. Tabl^Dlists key properties for three galaxies across the 
mass range that are analyzed in more detail; the complete list of 
NIHAO galaxies can be found in Wang et al. (2015). 

The hydrodynamical simulations are evolved using an im¬ 
proved version of the SPH code GASOLINE (Wadsley et al. 2004). 
The code includes a subgrid model for turbulent mixing of met¬ 
als and energy (Wadsley et al. 2008), heating and cooling include 
photoelectric heating of dust grains, ultraviolet (UV) heating and 
ionization and cooling due to hydrogen, helium and metals (Shen 
et al. 2010). 

For the NIHAO simulations we have used a revised treat¬ 
ment of hydrodynamics as described in Keller et al. (2014) that 
we refer to as ESF-GASOLINE2. Most important is the Ritchie & 
Thomas (2001) force expression that improves mixing and short¬ 
ens the destruction time for cold blobs (see Agertz et al. 2007). 
ESE-Gasoline2 also includes the time-step limiter suggested by 
Saitoh & Makino (2009), which is important in the presence of 


© 2015 RAS, MNRAS POQ.fTHTTI 





Core creation and destruction in galaxies 3 


strong shocks and temperature jumps. We also increased the num¬ 
ber of neighbor particles used in the calculation of the smoothed 
hydrodynamic properties from 32 to 50. 

2.1 Star formation and feedback 

The star formation and feedback modeling follows what was used 
in the MaGICC simulations (Stinson et al. 2013). Gas can form 
stars when it satisfies a temperature and a density threshold: T < 
15000 K and nth > 10.3 cm“®. Stars can feed energy back into the 
ISM via blast-wave supernova (SN) feedback (Stinson et al. 2006, 
see Agertz et al. 2013 for a discussion of the implementation of dif¬ 
ferent energy and momentum stellar feedback in galaxy formation 
simulations) and via ionizing radiation from massive stars before 
they turn in SN. Metals are produced by type II and type la SN. 
These, along with stellar winds from asymptotic giant branch stars 
also return mass to the ISM. The metals affect the cooling func¬ 
tion (Shen et al. 2010) and diffuse between gas particles (Wads- 
ley et al. 2008). The fraction of stellar mass that results in SN and 
winds is determined using the Chabrier (2003) stellar Initial Mass 
Function (IMF). 

There are two small changes from the MaGICC simulations. 
The change in number of neighbors and the new combination of 
softening length and particle mass means the threshold for star for¬ 
mation increased from 9.3 to 10.3 cm~^. The increased hydrody¬ 
namic mixing necessitated a small increase of pre-SN feedback ef¬ 
ficiency, eesf, from 0.1 to 0.13. This energy is ejected as thermal 
energy into the surrounding gas, which does not have its cooling 
disabled. Most of this energy is instantaneously radiated away, and 
the effective coupling is of the order of 1%. 

2.2 NIHAO galaxies stellar masses 

Since our aim is to study the impact of galaxy formation on the 
dark matter distribution it is very important to use realistic sim¬ 
ulated galaxies, i.e. galaxies able to reproduce the observed scal¬ 
ing relations. Haloes in our simulations were identified using halo 
finder AHeQ (Knollmann & Knebe 2009; Gill et al. 2004), adopting 
a spherical density contrast of 200 times the cosmic critical matter 
density. We dubbed the size of such a sphere i?vir, finally the stellar 
mass, Matar, is measured within a sphere of radius, Tgai = 0.2i?vir. 

Fig.E shows the stellar mass - halo mass relation for the 
NIHAO galaxies used in this paper, compared with the most 
commonly used abundance matching relations from the litera¬ 
ture: Brook et al. (2014), Kravtsov et al. (2014), Garrison-Kimmel 
et al. (2014), Behroozi et al. (2013) and Moster et al. (2013). This 
plot is similar to the one showed in Wang et al. (2015) with the only 
difference that here only the most resolved halo in each simulation 
is shown. As detailed in Wang et al. NIHAO galaxies are able to re¬ 
produce the the stellar mass - halo mass relation also at higher red- 
shift, and have realistic star formation rates for their stellar masses. 
Overall the unprecedented combination of high resolution and large 
statistical sample of the NIHAO suit offers a unique tool to study 
the response of DM to galaxy formation. 

2.3 Profile fitting 

To construct and fit the dark matter density profiles we used the 
same methodology as introduced in DC 14. Each Dark Matter den- 

^ http://popia.ft.uam.es/AMIGA 
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Figure 1. Stellar mass - halo mass relation for the NIHAO galaxies used in 
this work. All simulations have more than 400,000 particles in their virial 
radius (see Wang et al. 2015 for more details). The solid lines represent the 
most commonly used abundance matching results (see text). 


sity profile (in both hydro and N-body simulations) has been com¬ 
puted using regularly spaced shells on a logarithmic scale. The 
center of the halo has been determined using the shrinking sphere 
method (Power et al. 2003) and we set the minimum radius for the 
first shell as twice the softening length of the dark matter particles. 
Each density profile is computed up to the virial radius, which is 
determined using the spherical overdensity criterion with a virial 
overdensity of Avir = 200 times the critical density of the uni¬ 
verse. The error on the density in each bin is set according to the 
Poisson noise due to the finite number of particles each density bin. 

Following Di Cintio et al. (2014a) the DM central density 
slope (a) is subsequently fitted using a single power law, p cx r“, 
over a limited radial range, namely 0.01 < r/R^ir < 0.02, where 
Rvii is the virial radius. 


2.4 Profile convergence 

The determination of the minimum radius above which the results 
of a simulation are not affected by the finite resolution is a non 
trivial task (e.g. Gao et al. 2008). A quite commonnly used crite¬ 
rion for convergence has been suggested by Power et al. (2003) for 
collisionles simulations, and it is based on the two body relaxation 
time scale for particles in a gravitational potential. This empirical 
criterion radius ensures that the mean density inside the conver¬ 
gence radius is within 10 % of the value obtained in a simulation 
of much higher resolution (Schaller et al. 2015). 

Due to choice of having similar particle resolution in all the 
NIHAO galaxies, this convergence radius is quite constant across 
the whole mass range and at z = 0 is « 0.4 — 0.7% of Rvii- At 
z = 2, the inner region is marginally resolved if we use the origi¬ 
nal definition of Power et al. 2003, since we obtain a convergence 
radius of about 2% of the virial radius. On the other hand as no¬ 
ticed by Shaller et al. (2015), this definition is quite restrictive and 
it was based on pure DM simulations. For hydrodynamical sim¬ 
ulations it can be relaxed by requiring a 20% convergence in the 
enclosed density, (instead of 10%). If we follow the prescription of 
Schaller et al. the convergence radius is of the order 1.0-1.5% of 
Rvir, this means that while being “on the edge” we can still use the 
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r/Rvir 


Figure 2. Density profiles for three different simulations; on each plot we 
present the profile at 2 = 0 for the collisionless N-body (in black) and 
hydrodynamic (in red) simulation. The lines show power-law fits to the 
inner 1-2% of the halo, and are shown to 0. liJvir for clarity. The up¬ 
per panel shows an example where the DM halo contracts (halo mass of 
9.43 X IO^^Mq and stellar mass of 4.48 X IO^'^Mq), the middle panel 
shows a galaxy with a DM halo that expands (host halo mass of ), while 
the bottom panel shows a galaxy with very little change (halo mass of 
2.70 X and stellar mass of 4.28 X IO^Mq). 


same defintion for a up to 2 = 2, as long as we do not try to be too 
quantitative at these redshifts. 


3 RESULTS 

Thanks to the large range of masses covered hy our simulations we 
are able to see both expansion and contraction at 2 = 0 of the dark 
matter profile, depending on the halo mass and stellar content (see 
DC14 and references therein). 

An example of such a dichotomy is presented in Fig. [21 where 
we show a comparison of the dark matter density profile at 2 = 0 
in the hydro and in the N-body run for three galaxies with differ¬ 
ent total masses. The middle panel of Fig.|2]shows the DM density 
profile for a galaxy with a virial mass of « 2.2 x 10^^ M©. The pro¬ 
file in the hydro case (red points) shows a quite extended core. The 
situation is reversed for a more massive halo (« 9.4 x 
upper panel) where the DM profile is contracted in the hydro sim- 


Table 2. Best fit parameters for the value of a computed witihin 1 and 2 % 
of Rvir^S a function of Mstar/Mhalo- Afhalo and Mstar. 



n 

ni 

XQ 

Xi 

/3 

7 

M* 

Mh 

-0.158 

26.49 

8.77 xl0“® 

9.44x10“® 

0.85 

1.66 

Mh 

0.96 

7161 

9.96x10® 

1.85x10* 

1.01 

1.13 

M* 

0.53 

61.44 

2.48x10® 

6.83x10* 

29.5 

0.42 


ulation with respect to the N-hody run (hlack points). The lower 
panel shows a low mass galaxy (halo mass « 2.7 x IO^^Mq), in 
this case the both profiles are cuspy, but it is clear that the hydro 
simulation has a shallower central slope than the N-body one. 

This plot shows that simulations that include a realistic treat¬ 
ment of hydrodynamics do not share a universal dark matter den¬ 
sity profile in contrast with dark matter only simulations. We go 
into detail about how the universality is broken in section 4. 

As already pointed out by Di Cintio et al. (2014a) there is tight 
relation between the DM response and the star formation efficiency 
of a galaxy, defined as the ratio between stellar mass and halo mass. 
Fig. [3] shows the relation between a and the Mstar/Mhaio ratio. 
The red points with errorbars are the NIFIAO simulations results at 
2 = 0. The black open circles are the results from the N-body runs, 
and for these simulations we use the stellar mass of their hydrody- 
namical counterpart. 

The NIHAO results show a quite interesting behaviour: the 
inner slope, a, varies as a function of star formation efficiency. 
At low values of Mstar/Muaio the N-body and hydro results are 
quite similar and both predict a constant value of a ~ —1.5. Then 
a increases steadily and reaches a maximum for Mstar/Mhaio ~ 
6 X 10“®, in agreement with previous findings from DC14 (black 
dashed line). At the highest masses, after the maximum expansion, 
a decreases and for very large values of Mstar/Mhaio of the order 
of few percent, the hydro slopes become steeper than their N-body 
counterparts, suggesting halo contraction on those scales. This is in 
agreement with previous results from Di Cintio et al. 2014b, where 
it was also found that at the mass scales of the Milky Way, the 
DM profiles have a much higher concentration parameter in hy¬ 
dro simulations compared to collisionless results. It is also worth 
remembering that, as shown in Fig. [T] all our galaxies are on the 
observed stellar mass-halo mass relation (e.g. Behroozi et al. 2013; 
Kravtsov et al. 2014), and hence our halo contraction is not due 
to unphysical overcooling as it happened in older simulations (e.g. 
Gnedin et al. 2004) and it is consistent with recent results from 
other groups on the same mass scale (Schaller et al. 2015). 

We tried to capture the behaviour of the slope of the den¬ 
sity profile as a function of the star formation efficiency (x = 
Mstar/Mhaio) with a different fitting formula than DC14: 


a{x) 


n - logio 


ni 




( 1 ) 


This new fitting function has the advantage that it converges 
to a fixed value (« —1.5) for x that goes to zero, in contrast to the 
power-law behavior of the initial suggestion from DC14. Since, for 
continuity, we do expect to recover the slope of the N-body results 
for a negligible amount of stars {x < 10“®), this new formula 
provides a better fit to the simulation data. The grey area is the one 
sigma scatter around the mean value of a of « 0.18. 
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Figure 3. Inner dark matter density slope, a, as a function of Mstar /Afjjaio for galaxies at z = 0. The filled red circles show results for hydro runs (with 
errorbars from the profile fit) while the black open circles show results for the N-hody DM only runs. The dotted line is the original function from Di Cintio 
et al. (2014a), the thick line is the refitted curve with our new simulations using the Planck cosmology with the shaded region showing the la scatter of 0.18. 


The difference between the two curves for Mstar/Afhaio > 
10“® is mainly due to the different cosmological model adopted 
in the two studies (Planck vs. WMAP3), while other differences, 
as for example the improved hydrodynamics, are just a secondary 
effects (see Stinson et al. 2015 for a more detailed comparison be¬ 
tween MaGICC and NIHAO galaxies). In the Planck cosmology 
the halo concentrations are higher by « 45% (Dutton & Maccio 
2014), which explains why we get lower (cuspier) a-values with 
respect to DC14. Table[3]lists the values of our new relation for a 
vs. Mstar/Mhaio, (updated to the Planck Cosmology). 

In Fig.|4]we show the same relation as in Fig.j^but at redshift 
2 = 1 (upper panel) and at redshift 2 = 2 (lower panel). The line 
in the plot is still the 2 = 0 fitting result. As the galaxy evolve, dark 
matter accretion and star formation generally move points from left 
to right. The points still cluster around the 2 = 0 relation, but scat¬ 
ter increases with redshift, so integrated star formation efficiency 
becomes a less useful metric when the galaxies are young. 

The origin of the scatter is two fold. Part of it is due to resolu¬ 
tion effects, at higher redshift the virial radius is smaller, this means 
that we are probing regions closer to our resolution limit which im¬ 
plies a lower number of particles in the inner region of the halo 
and hence a larger Poisson noise. Part of the scatter is also physical 
and related to the bursty nature of the star formation history at high 
redshift, that induces short time variations in the density profile (see 
next section). 

To check whether star formation efficiency best correlates 
with a in our large sample of galaxies. Fig. |5] and Fig. |6] show 
how a varies with the halo mass and the stellar mass of our galax¬ 
ies respectively. In both cases the behaviour of a is similar to the 
one seen in Fig.O with a maximum in ’core’ creation around 10*^^ 


(10®) Mq in halo (stellar) mass. These plots show a similar scat¬ 
ter in Q as in Fig. for low values of the halo (stellar) mass: 
Mhaio < 5 X (Mstar < lO^M©), but a larger scatter 

for a at high masses. This confirms the earlier results of Di Cin¬ 
tio et al. (2014a) where using galaxies run with different stellar 
feedback prescription, they show that the integrated star formation 
efficiency is the best parameter to capture the effect of baryons on 
the DM distrubution. 

The role of the integrated star formation efficiency is also im¬ 
portant to explain the similiraties and (small, partial) differences 
between our results and some previous studies on the same topic 
from different groups. 

Madau et al. (2014) analyzed the evolution of the density 
profile of a small group of dwarf galaxies, with only 4 of them 
with containing stars. The simulations were run with an earlier 
version of the GASOLINE code, with similar SN feedback. Two 
of those galaxies (dubbed Doc and Bashful) have masses around 
1 — 3 X 10^°Mq and do show extended cores, somehow in con¬ 
trast with the results shown in figure |5] On the other hand the 
galaxies presented by Madau and collaborators have over-massive 
stellar bodies, that over predict results from Local Group abun¬ 
dance matching from Brook et al. (2014) and Garrison-Kimmel 
et al. (2014) by a factor > 10 for “Doc” and of about 3 — 4 
for “Bashful” (w.r.t Garrison-Kimmel et al. even more for Brook 
et al.). This overcooling can be explained by the lack in their simu¬ 
lation of any other source of feedback besides SuperNovae (see for 
example the results of Stinson et al. 2013). When the higher (may 
be too high) stellar mass -halo mass ratio of these galaxies, com¬ 
pared to our NIHAO sample, is taken into account the results are in 
very good in agreement. 
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Figure 4. Inner dark matter density slope, a, as a function of Mstar/ Afhaio 
for galaxies at 2 = 1 (upper panel) and at 2 = 2 (lower panel). The solid 
line and shaded region is 2 : = 0 relation from Fig. [5] 



Figure 5. Inner dark matter density slope, a, as a function of the halo mass 
for galaxies at 2 : = 0. The solid line and shaded region are the fitted relation 
using Eq.^and its scatter. The parameters of the fitting are reported in Table 

m 

A similar agreement is also found when comparing our rusults 
with recent simulations presented by the FIRE collaboration in 
Chan et al. (2015). For example there is a quite remarkable sim¬ 
ilarity in the peak of the core formation and in the non monotonic 
response of DM to galaxy formation. This similarity is even more 
striking when we take into account that the two simulation sets (NI- 
HAO and FIRE) differ in resolution, numerical technique and stel¬ 
lar feedback implementation. 

One place where our and their results do not agree is in the 
decline of the value of a at low masses and which brought them 
to suggest the possibility of a “threshold” halo mass at around 
10^° Mq needed to develop large cores. It is worth mentioning that 
while the EIRE simulations reach a much higher resolution than our 
current sample, they don’t have large number statistics, with only 9 
galaxies across 5 orders of magnitude in stellar mass. 



Figure 6. Inner dark matter density slope, a, as a function of the stellar 
mass for galaxies at z = 0. The solid line and shaded region are the fitted 
relation using Eq.[T]and its scatter. The parameters of the fitting are reported 
in Table |3] 

We do not see such a threshold in our sample at this halo mass 
(but note that we do see a substantial change in a for a stellar mass 
of « Mq ). This difference can have many explanations. The 
higher resolution of the EIRE simulations might model the inter¬ 
play between star formation and the outflows they drive better The 
partial evidence for a threshold could be due to the much lower 
sample size of Chan et al. : they only have three galaxies in the 
dark matter range between 5 x 10® and 5 x 10 ^®Mq, while we 
have twenty five, and hence we are less subject to artificial “jumps” 
due to low number statistics and intrinsic scatter. 

Einally as in the case of the Madau et al. results, the FIRE sim¬ 
ulations predict a stronger star formation efficiency at low masses, 
which partially over produces abundance matching results from 
Garrison-Kimmel et al. (2014) and Brook et al. (2014). This de¬ 
parture from abundance matching results, starts exactly above halo 
masses of 10^® Mq (Dutton et al. 2015, figure 1). Overall we would 
like to stress again the importance of having large samples and re¬ 
alistic stellar to halo masses ratios in order to infer the behavior of 
the dark matter distribution in hydrodynamic simulations. 

3.1 Single profile evolution 

To study the evolution of individual halos, we make case studies of 
3 simulations, g2.63el0, g2.19ell, and g8.06ell (already shown 
in Eig.l^and presented in table[T), that show different behaviors as 
a function of redshift. In Fig.|7] we show the DM density profiles at 
2 = 0 (red points) and 2 = 2 (green points) for the three galaxies. 
At 2 = 2 the total halo masses are 1.48 x 2.61 x 10^® and 
1.50 X 10*^^, while the stellar masses are 3.29 x 10^, 2.64 x lO’^, 
and 1.11 X 10®, for g2.63el0 , g2.19el 1, and g8.06el 1 respectively. 
The dashed line shows a power-law fit to the density between 1 and 
2% of the virial radius. 

As expected, at 2 = 2 the middle mass galaxy already has 
an expanded profile (compared to the N-body run), with a slope of 
a = —0.65, and then the inner density flattens by the end of the 
simulation. At 2 = 2, the dark matter profile of the higher mass 
halo looks nearly identical to the middle mass halo at 2 = 0. It 
evolves differently; its central density profile becomes steeper with 
time. The low mass galaxy (g2.63el0) has a very mild evolution 
and the slope of the profiles doesn’t change much from 2 = 2 to 
the present day. 

This strong time variation of the size of the central den- 
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r/Rvir 


Figure 7. DM density profiles for our three test galaxies at two redshifts. 
The profiles at 2 = 0 are represented in red, and at 2 = 2 in green. The 
dotted line is fitted for each profile between one and two percent of the virial 
radius and is shown up to ten percent of the virial radius to help visualize. 


sity core are correlated to quick variations in the Star Formation 
Rate of the galaxy (Pontzen & Govemato 2012, see also Maccio 
et al. (2012) for the effect of the sloshing of the gas inside the po¬ 
tential well). 

Fig- m provides example cases of how bursty star formation 
effects the inner dark matter profile slope for galaxies of three dif¬ 
ferent masses spanning the NIHAO sample. The lowest mass case, 
g2.63el0, has early bursts of star formation before its star formation 
dies away entirely. The bursts of star formation are not reflected in 
the inner profile slope at early times, most likely because of the low 
resolution of the profiles that makes it hard to measure the slope. 
As time goes on and the resolution increases, the slope stays mod¬ 
erately flatter than NFW, at a value slightly higher than —1. 

The middle mass case, g2.19ell, displays bursts of star for¬ 
mation throughout its evolution. The largest bursts start around 7 
Gyr {z ~ 1). The bursts coincide with a flattening of the inner dark 
matter profile slope. 

The highest mass galaxy, gS.Ohell, shows a monotonically 
rising star formation history. Throughout the history, there are 
bursty periods. Its inner density profile evolves from a cusp to a 
core before contracting back into a cusp with a. = —1.7 at 2 = 0. 



1 2 3 4 5 6 7 8 9 10 11 12 13 


Time [Gyr] 

Figure 8. Time evolution of the Star Formation Rate (SFR, upper panels) 
and the central density slope (a, lower panel), for our three test galaxies: 
g8.06ell (blue), g2.19ell (red), g2.63el0 (orange). 

It changes from a cusp to a core at the same time that its star for¬ 
mation rate increases by a factor of 4, around 7 Gyr. 

3.2 Core creation 

The variation of the central slope of the density profile is also linked 
to the creation (and destruction) of a constant density inner core. In 
order to have an estimate of the size of such a core, we have decided 
to follow the approach of Maccio et al. (2012) and fit the profiles 
with the following parametric description, originally presented in 
Stadel et al. (2009): 

p(r) =poexp(A[ln(l-|-r/i7A)]^). (2) 
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4 5 6 7 8 9 10 11 1213 



Figure 9. Core size Rx as a function of time 1/(1 + z) for g2.19el 1 (red) 
and g8.06ell (blue). The grey shaded region corresponds to the resolution 
limit of the simulations. 

In this parameterization, the density profile is linear down to a 
scale R\ beyond which it approaches the central maximum den¬ 
sity po as r —>■ 0. We also note that if one makes a plot of 
dlnp/dln(l -f r/Ry) versus ln(l -f r/Ry) then this profile forms 
an exact straight line with slope 2A. This fitting function is ex¬ 
tremely flexible and makes it possible to reproduce at the same time 
both cuspy and cored density profiles. In the following we will use 
the fitted value of R\ as the core size. 

In Fig.l^we show the core size (R\) as function of time for the 
two more massive galaxies shown in Fig. [7] The grey band in the 
lower panel of the plot represents the softening of the simulations, 
any core value below this limit has to be interpreted as a non-cored 
profile. The lower mass galaxy (red line) shows a gradually increas¬ 
ing average core size from 2 = 2 to 2 = 0. The behavior of the 
more massive galaxy (blue line) is much more interesting. At high 
redshift (2 = 2) the halo has a core of order ~ 1.5 kpc, in between 
2 = 1.5 and 2 = 1 the core size fluctuates a lot and reaches quite 
large value up to 4 kpc. The core is then definitely destroyed after 
2 = 0.5 and the galaxy shows no presence of a central constant 
density, as confirmed by the profile shown in Fig.[2 

It is interesting to look deeper into the origin of the spikes of 
the core value around 2 ~ 0.5 (T ~ 8.5 Gyrs) for the g8.Obeli 
halo. In Fig. [H we show the dark matter density profile of this 
galaxy at the time of the maximum size of the core (2 = 0.67, 
black line), and right before the core creation (red line) and right 
after the core destruction (blue line). It is clear from the figure that 
the profile does change quite strongly in the time interval spanned 
by the spike. A large core (up to 4 kpc) is created and then destroyed 
in less than a Gyr. 

This fluctuation are related to the rapid increas of an order of 
magnitude of star formation from T=6.5 Gyrs to T=8.5. The mor¬ 
phological analysis of the galaxy shows that this extended burst is 
due to disc instability, possibly triggered by a satellite accretion. 
During this phase a lot of gas is both accreated and ejected from 
the central region creating the perfect conditions for core creation 
and destruction (e.g. a la Pontzen & Govemato 2012). 


3.3 Core destruction 

Our simulations clearly show that dark matter halo profiles are not 
a static entity. On the contrary they change significantly during the 
formation and evolution of the halo. While all haloes start as being 
cuspy, as predicted by pure CDM simulations, they can develop a 



Figure 10. Evolution of the Dark Matter density profile for the g8.Obeli 
halo around 2 ~ 0.7. The profile changes significantly with clear evidence 
for the creation and subsequent destruction of a large density core. The time 
spanned by the three profiles is of the order of one Gyr. The two dotted lines 
mark the region within which a is computed. 


core. Depending on the stellar to halo mass ratio, the core can per¬ 
sist to the present day, or be destroyed, being replaced by a central 
cusp. 

In massive haloes (e.g. with halo mass around lO^^M©) the 
cusp regeneration is mainly due to strong gas inflow in the central 
region, which builds up a large stellar body in the center of the 
galaxy causing a deepening of the local potential well. We will now 
analyze these processes in more detail. 

In Fig.[TT]we show the relative abundance of dark matter, gas 
and stars in the inner 2 kpc for our test galaxies. In the case of 
the lowest mass galaxy (g2.63el0) the central potential is always 
dominated by dark matter which accounts for more than ~ 90% 
of the total mass in the center. The lack of evolution in a shown 
in Fig. [8] is then due to the very “passive” evolution of the central 
region of this halo after star formation is over. Nevertheless the final 
effect of baryons on this low mass halo is still a net expansion (as 
shown in figure^ with a = —0.9 in the Hydro run with respect to 
a value of ~ —1.6 in the N-body one; despite stars and gas account 
for only ~ 10% of the mass inside 2 kpc. 

For the middle mass galaxy (which retains a cored profile all 
the way to 2 = 0), the potential in the central region is dominated 
by dark matter and gas, while stars are sub-dominant. The amount 
of gas is rapidly varying with time and this causes quick potential 
variations and a rapid movement of the center of mass of the gas 
compared to that of the dark matter; both these effects contribute in 
creating and maintaining a density core. 

The situation is quite different for the most massive galaxy: 
after T = 8 Gyr the potential in central region is dominated by 
stars with the gas (and DM) being highly sub-dominant. At this 
point, the deeper potential well caused by the stars is anchoring the 
global potential and making any possible gas outflow a negligible 
effect. The dark matter reacts to the new, enhanced central potential 
by contracting towards the center. This contraction is causing the 
core size to drop below our resolution (Fig. H and the profile slope 
a to drop below -1.5 (Fig. [8j. 

The deepening of the central potential (due to gas inflow and 
star formation) in our more massive halo can be seen in Fig. M 
where we plot the evolution of the total mass (DM-l-gas-l-stars) 
within different radii for our three galaxies. In the middle mass 
halo (g2.19ell) the mass is constant practically at all radii after 
8 Gyrs, this implies a constant gravitational potential. In the higher 
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Time [Gyr] 


Figure 11. The contribution of DM (black), gas (blue) and stars (red) to the 
total mass in the inner 2 kpc for g8.Obeli (upper panel) g2.19ell (middle 
panel) and g2.63el0 (lower panel), as a function of time. Note the log y-axis 
in the lower panel. 

mass halo (g8.06ell) at large radii (100 kpc, still well within the 
virial radius) the total mass is also constant, but at small radii (5 
and 10 kpc) the mass rapidly increases by a factor of three after t=8 
Gyrs. This higher mass caused the gravitational potential to deepen 
and reduces the ability of baryonic feedback to counteract the DM 
contraction. 

It is worth noticing that the change in the mass distribution is 
mainly local (i.e. limited to the regiond dominates by baryons) and 
not global, something that should be taken into account when com¬ 
puting analytic estimates of the total energy input from SN needed 
to alter the DM distribution. This is a still quite debated topic in the 
literature and our simulations can be used to determine the valid¬ 
ity of different analytic approaches (e.g. Penarrubia et al. 2012 and 
Maxwell et al. 2015). 


4 COMPARISON WITH OBSERVATIONS 

There is a very large amount of literature on the DM density pro¬ 
file in observed galaxies (e.g Moore 1994, Salucci & Burkert 2000; 
Dutton et al. 2005; Simon et al. 2005; de Blok et al. 2008; Kuzio 
de Naray, McGaugh & de Blok 2008; Kuzio de Naray, McGaugh 



2 4 6 8 10 12 14 

Time [Gyr] 


Figure 12. Time evolution of the tola/ mass enclosed within different radii, 
namely 5,10,50 and 100 kpc, normalized to the z = 0 mass. From top 
to bottom the results for our three test galaxies: g8.Obeli, g2.19ell, and 
g2.b3el0. The strong evolution of the inner mass (2 and 5 kpc) for galaxy 
g8.Obeli is thought to be the cause of the halo recontraction (see text for 
more details). 

& Mihos 2009; Oh et al. 2011a, Oh et al. 2015). Despite the quite 
large scatter in the slope of the observed profiles there is a clear 
indication for aobs > —1 (e.g. Oh et al. 2011a), which is in con¬ 
trast with naive expectations from collisionless simulations. This 
disagreement is strongly alleviated when the effects of baryons 
are taken into account, as recently shown in Karukes, Salucci & 
Gentile (2015), where the authors compared the measured DM 
density profile in NGC 3198 with the predictions from the DC14 
model. Thanks to the large number of simulated galaxies in the NI- 
HAO project we can for the first time extend previous observation- 
simulation comparisons that were based on a limited number of 
objects (e.g. Oh et al. 2011b). 

In Fig. [T^ we show the inner slope of the dark matter density 
profile a vs. the radius Rin of the innermost point within which 
a is measured for both observations and simulations. The obser¬ 
vational data are split in three groups, the grey symbols are results 
from dwarf and Low Surface Brightness (LSB) galaxies (de Blok et 
al. 2001; de Blok & Bosma 2002; Swaters et al. 2003), cyan filled 
squares are the results from the THINGS survey (Oh et al. 201 la), 
while open blue squares with error-bars are the results from the 
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Figure 13. The inner slope of the dark matter density profile a vs. the ra¬ 
dius iJin of the innermost point within which a is measured. Grey symbols 
are resulls from LSB galaxies (de Blok el al. 2001; de Blok & Bosma 2002; 
Swaters et al. 2003), cyan filled squares are the results from the THINGS 
survey (Oh et al. 2011a), open blue squares with error-bars are the results 
from the LITTLE THINGS survey (Oh et al. 2015). The red points are all 
NIHAO galaxies wilh stellar masses between 10^ — to mimic the 

range of the observational points. The Isothermal profile line is for a core 
of 1 kpc, while the NEW and Einasto lines are for a concentration of 14 
(Dutton & Maccio 2014). 


LITTLE THINGS survey (Oh et al. 2015). The three lines show 
theoretical predictions from a given profile: the dotted line is for 
an isothermal profile with a core size of 1 kpc, the solid line is for 
the Einasto profile (Einasto 1965, Merritt et al. 2005) with a con¬ 
centration of 14 as expected for dark matter haloes around 10^^ Mq 
(Dutton & Maccio 2014) and the dashed line is for the NEW profile 
with the same concentration as the Einasto one. It is interesting to 
note that on the radial scales of this plot, the Einasto profile (which 
is known to provide a better fit to simulations, Merritt et al. 2005, 
Navarro et al. 2010) predicts a steeper slope for DM haloes with 
respect to NEW. 

The red circles with error-bars are the predictions for a from 
the NIHAO simulation suite at 1-2% of Rvii, for galaxies with stel¬ 
lar masses between 10^ — IO^^Mq to be consistent with the LIT¬ 
TLE THINGS sample range (Oh et al. 2015). All the NIHAO dark 
mater haloes are expanded with respect to collisionless simulations 
expectations. We don’t have enough resolution to probe the inner 
part of the radial profile at the same extent as the LITTLE THING 
survey, nevertheless our results are in very good agreement with 
the observations where they overlap in Rin- Moreover thanks to the 
large number of simulated galaxies, we are also able to reproduce 
the observed scatter in q at a fixed radius. 

This plot further demonstrates the importance of considering 
the effect of baryons on the dark matter distribution when compar¬ 
ing observations and simulations based on the Cold Dark Matter 
model. 


5 CONCLUSIONS 

We have used the NIHAO simulation suite (Wang et al. 2015) to 
study the effects of baryons on the radial distribution of dark matter 
in collapsed haloes. In comparison with a similar study performed 
earlier with the MaGICC sample (Di Cintio et al. 2014a, DC 14) we 
have a larger and more coherent set of simulations, which have been 
performed with the latest compilation of cosmological parameters 
from the Planck Collaboration (2014). Moreover we use a substan¬ 
tially improved version of the GASOLINE code, which thanks to a 
new formulation of the SPH force computation (Keller et al. 2014), 
aims to solve the known short comings of the classical SPH imple¬ 
mentation (Agertz et al. 2007; Hopkins 2013). 

Overall we confirm previous findings from DC 14, namely that 
the effect of baryons on the slope (a) of the central distribution of 
dark matter strongly depends on the ratio between the stellar mass 
and total mass of the halo. At the present time, haloes with an in¬ 
termediate ratio of stellar-to-halo mass (10~® < Mstar/AThalo < 
10“^ tend to expand their profiles with respect to pure collisional 
simulations, haloes with lower values of Mstar/Afhaio tend to re¬ 
tain the original dark matter cuspy profile, while haloes with large 
(ATstar/AThaio > 10“^) ratio exhibit contracted profiles. There is 
a smooth transition among these different behaviors and we have 
parameterize it with a new, simple but accurate fitting formula. The 
slopes predicted by the NIHAO simulation suite are in excellent 
agrement with high resolution observations of dwarf galaxies from 
the THINGS and LITTLE THINGS surveys (Oh et al. 2015). 

When the analysis is, for the first time, extended to higher red- 
shifts, we find that haloes tend to cluster around the redshift zero a 
vs. Mstar/AThaio relation, even though the scatter around the mean 
increases with redshift. Our findings imply thal even haloes that 
have very cuspy profiles loday (as for example for Milky Way mass 
haloes), likely had a cored dark matter distribution in the past. 

We carefully looked at three specific haloes with masses 
around 10^°, 10^^, lO^^M©. The lightest shows a central slope of 
about —0.9 at 2 = 0, which is consistent with a cuspy profile, but 
is also very much shallower than what a pure N-body simulations 
predicts at the same mass and radii {a ^ —1.6), pointing to a clear 
halo expansion even at this low mass scales. The density profile 
was even shallower at higher redshift {z ~ 1) when the halo went 
through several star formation episodes. 

The other two haloes present a quite dynamic evolution in 
their central regions as a function of time, showing that not only 
cores are created and destroyed during the galaxy evolution but that 
this phenomenon happens on relatively (cosmologically) short time 
scales of the order of a Gyr. 

Overall our results show that one of the most firm predictions 
of the Cold Dark Matter theory, namely the existence of a “univer¬ 
sal profile”, is no longer valid when galaxy formation is taken into 
account. The central part of any dark matter halo has had a much 
more tumultuous life than what dissipationless simulations predict. 
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